      SUBROUTINE  ROBSWF(CH,XI,M,N,MINT,INDX,NMAX,MAXJW,JW,D,ALAMD,     00010000
     +                   SIGMA,DSUM,R1,DR1,R2,DR2)                      00020011
CCC                                                                     00030000
CCC       THE OBLATE SPHEROIDAL RADIAL FUNCTIONS WITH REAL ARGUMENTS.   00040000
CCC       THE SIZE PARAMETER 'CH' IS REAL.                              00050000
CCC       THE FIRST KIND FUNCTIONS AS A SERIES OF THE SPHERICAL BESSEL F00060000
CCC       INDX=1, EXTERNAL FIELD OR C1 ;  INDX=2, INTERNAL FIELD OR C2. 00070000
CCC       MINT=INDEX PARAMETER FOR THE CALCULATION OF BESSEL FUNCT.     00080000
CCC       THE SECOND KIND FUNCTIONS ARE NEVER CALCULATED FOR INDX = 2.  00090000
CCC                                                                     00100000
      IMPLICIT REAL*8(A-H,O-Z)                                          00110000
      PARAMETER  (JMX=150,MXPQ=190,MXSJ=150,MXSN=190)                   00120006
      COMPLEX*16 CI,C,C2,XX,QNA,QNB,DQN,QD,QSUM,DQSUM,PSUM,DPSUM,COEF2, 00130000
     1           S1,S2,S3,B,HSUM,DHSUM,HCONV,R3,DR3,QN(80),             00140008
     2           PX(MXPQ),DPX(MXPQ),QX(MXPQ),DQX(MXPQ)                  00150002
      DIMENSION D(JMX),DN(80),DP(200),DNRAT(80),DPRAT(200),AMX(JMX),    00160006
     1          SJ(MXSJ,2),DSJ(MXSJ,2),SN(MXSN,2),DSN(MXSN,2),          00170002
     2          UP(MXPQ),DW(500),A(MXPQ),CR(3),R2M(3),DR2M(3)           00180006
      FLOAT(N)=DFLOAT(N)                                                00190000
      E(J)=C2*DFLOAT((M+M+J)*(M+M+J-1))/DFLOAT((2*(M+J)-1)*(2*(M+J)+1)) 00200000
      F(J)=DFLOAT((M+J)*(M+J+1))+C2*DFLOAT(2*(M+J)*(M+J+1)-2*M*M-1)/    00210000
     +     DFLOAT((2*(M+J)-1)*(2*(M+J)+3))-ALAMD                        00220002
      G(J)=C2*DFLOAT((J+1)*(J+2))/DFLOAT((2*(M+J)+1)*(2*(M+J)+3))       00230000
CC    
C     CALL ERRSET(207,256,-1,0,0,208)                                   00240024
CC                                                                      00250000
      PAI2=3.141592653589793D0/2.D0                                     00260002
      CI=(0.D0,1.D0)                                                    00270000
      C=(0.D0,-1.D0)*CH                                                 00280000
      C2= -CH*CH                                                        00290000
      CH2= CH*CH                                                        00300002
      XI2=XI*XI+1.D0                                                    00310000
      XX=CI*XI                                                          00320000
      H=CH                                                              00330000
      HXX=H*XI*XI                                                       00340000
      L=N-M                                                             00350000
C                                                                       00360014
        IF(M.EQ.0)  THEN                                                00370002
      XIA=1.D0                                                          00380002
      XIB=0.D0                                                          00390002
        ELSE                                                            00400002
      XIA=DSQRT(XI2/(XI*XI))**M                                         00410000
      XIB=-DFLOAT(M)/(XI*XI2)                                           00420000
        END IF                                                          00430002
      IF(MAXJW.GE.JMX-1)  MAXJW= JMX-2                                  00440007
CC             CRITERION FOR THE METHOD OF COMPUTATION                  00450002
      IF(XI.GE.0.80) GO TO 5000                                         00460000
      ALP=0.00667*H+0.7334                                              00470000
      LM0=0.4*H-0.999                                                   00480000
      IF(XI.GE.0.5) LM0=0.6*H-1.999                                     00490000
      GO TO 5010                                                        00500000
 5000 LM0=H-0.999                                                       00510000
      IF(H.GT.25.0) LM0=1.9*H-23.499                                    00520000
      ALP=0.00667*H+0.73334                                             00530000
      IF(H.GT.28.0) ALP=1.3                                             00540000
 5010 MTHD=LM0-ALP*FLOAT(M)                                             00550000
      IF(MTHD.LT.0)  MTHD=0                                             00560000
      LM1=MTHD-4                                                        00570000
      LM2=MTHD+4                                                        00580000
      LMTN=LM2+4                                                        00590000
      IF (M.GE.20.AND.H.GE.20.0) LMTN=LM2+FLOAT(M)*H/25.                00600000
      CR(1)= 1.D+25                                                     00610000
      CR(2)= 1.D+10                                                     00620000
      CR(3)= 1.D+15                                                     00630000
      R2=0.D0                                                           00640000
      DR2=0.D0                                                          00650000
C                                                                       00660002
C         SPHERICAL BESSEL FUNCTIONS                                    00670000
C                                                                       00680002
      IF(M.EQ.MINT.AND.N.EQ.M) GO TO 4002                               00690000
      GO TO 4005                                                        00700000
 4002 NSJ=MAXJW+M                                                       00710000
        IF(NSJ.GE.MXSJ-1)  NSJ= MXSJ-2                                  00720009
      DX=CH*XI                                                          00730000
      W=( DSIN(DX)-DX* DCOS(DX))/(DX*DX)                                00740000
      T3= 0.D0                                                          00750000
      T2= 1.D-70                                                        00760000
      NUP=10+DX/5.                                                      00770000
      NSWT=NSJ                                                          00780000
      NA=NSJ+NUP                                                        00790000
      DO 4000 NB=1,NA                                                   00800000
      NC=NA+1-NB                                                        00810000
      T1=DFLOAT(NC+NC+1)*T2/DX-T3                                       00820000
      IF( DABS(T1)-1.D70) 4001,4003,4003                                00830000
 4003 T1=T1*1.0D-70                                                     00840000
      T2=T2*1.0D-70                                                     00850000
      NSWT=NC                                                           00860000
 4001 T3=T2                                                             00870000
      T2=T1                                                             00880000
      IF(NC.GT.NSJ) GO TO 4000                                          00890000
      SJ(NC,INDX)=T1                                                    00900000
 4000 CONTINUE                                                          00910000
      RATIO=W/SJ(2,INDX)                                                00920000
      DO 4010 ND=1,NSJ                                                  00930000
      IF(ND.GT.NSWT) GO TO 4015                                         00940000
      SJ(ND,INDX)=RATIO*SJ(ND,INDX)                                     00950000
      GO TO 4010                                                        00960000
 4015 SJFACT=1.D-70                                                     00970000
      IF( DABS(SJ(ND,INDX)).LT.1.0D-5) SJFACT=0.D0                      00980000
      SJ(ND,INDX)=RATIO*SJ(ND,INDX)*SJFACT                              00990000
 4010 CONTINUE                                                          01000000
      DSJ(1,INDX)=-SJ(2,INDX)                                           01010000
      DO 4020 NE=2,NSJ                                                  01020000
      NF=NE-1                                                           01030000
 4020 DSJ(NF+1,INDX)=SJ(NF,INDX)-DFLOAT(NF+1)*SJ(NF+1,INDX)/DX          01040000
C                                                                       01050002
C         SPHERICAL NEUMANN FUNCTIONS                                   01060002
C                                                                       01070002
      IF(INDX.EQ.2) GO TO 4005                                          01080000
      IF(XI.LT.1.005) GO TO 4005                                        01090000
      IF(HXX.LT.15.0) GO TO 4005                                        01100000
      SCLSN=1.0E-20                                                     01110000
      IF( DABS(DX).GE.15.0) SCLSN=1.0E-10                               01120000
      SN(1,INDX)=- DCOS(DX)/DX                                          01130000
      SN(2,INDX)=-( DCOS(DX)+DX* DSIN(DX))/(DX*DX)                      01140000
      SN(1,INDX)=SN(1,INDX)*SCLSN                                       01150000
      SN(2,INDX)=SN(2,INDX)*SCLSN                                       01160000
      DSN(1,INDX)=-SN(2,INDX)                                           01170000
      MAXSN=MAXJW+M                                                     01180000
       IF(MAXSN.GE.MXSN-1)  MAXSN= MXSN-2                               01190006
      DO 4030 NP=2,MAXSN                                                01200000
      NQ=NP-1                                                           01210000
      IF(NQ.LE.1) GO TO 4040                                            01220000
      SN(NQ+1,INDX)=DFLOAT(NQ+NQ-1)*SN(NQ,INDX)/DX-SN(NQ-1,INDX)        01230000
 4040 DSN(NQ+1,INDX)=SN(NQ,INDX)-DFLOAT(NQ+1)*SN(NQ+1,INDX)/DX          01240000
 4030 CONTINUE                                                          01250000
 4005 CONTINUE                                                          01260000
C                                                                       01270002
C         ASSOCIATED LEGNDRE FUNCTIONS                                  01280000
C                                                                       01290002
      IF(INDX.EQ.2) GO TO 4007                                          01300000
      IF(N.GT.M) GO TO 4007                                             01310000
      MAXPQ=MAXJW+M                                                     01320000
       IF(MAXPQ.GE.MXPQ-1)  MAXPQ= MXPQ-2                               01330006
      MLEG=0                                                            01340000
      ALEG=1.D0                                                         01350000
 4070 MLEG=MLEG+1                                                       01360000
      ALEG=ALEG*DFLOAT(2*MLEG-1)                                        01370000
      IF(MLEG.LT.M)  GO TO 4070                                         01380000
         CALL  CLEGQ(XX,M,MAXPQ,QX)                                     01390002
      PX(1)=ALEG*(CDSQRT(XX*XX-1.D0)**M)                                01400000
      PX(2)=DFLOAT(2*M+1)*XX*PX(1)                                      01410000
      DO 4050 I=1,MAXPQ                                                 01420000
      MI=M+I                                                            01430000
      IN=I-1                                                            01440000
      PX(I+2)=(DFLOAT(2*MI+1)*XX*PX(I+1)-DFLOAT(M+MI)*PX(I))/DFLOAT(I+1)01450000
      DPX(I)=(DFLOAT(I)*PX(I+1)-DFLOAT(M+I)*XX*PX(I))/(XX*XX-1.0)       01460000
      DPX(I)=CI*DPX(I)                                                  01470000
      IF(I.EQ.1) GO TO 4050                                             01480000
      DQX(I-1)=(DFLOAT(IN-M)*QX(I)-DFLOAT(IN)*XX*QX(I-1))/(XX*XX-1.0)   01490000
      DQX(I-1)=CI*DQX(I-1)                                              01500000
 4050 CONTINUE                                                          01510000
          WRITE(6,888)  (KP, PX(KP),KP=1, MAXPQ)                        01511024
  888     FORMAT(1H / 5(I4,1P2E11.3) )                                  01512027
      IF(M.EQ.0) GO TO 4007                                             01520000
      S1=QX(2)                                                          01530000
      S2=QX(1)                                                          01540000
      DO 4080 KA=1,M                                                    01550000
      KK=KA-1                                                           01560000
      S3=(DFLOAT(1-2*KK)*XX*S2+DFLOAT(KK+M-1)*S1)/DFLOAT(M-KK)          01570000
      S1=S2                                                             01580000
      S2=S3                                                             01590000
      QN(KA)=S3                                                         01600000
 4080 CONTINUE                                                          01610000
 4007 CONTINUE                                                          01620000
C                                                                       01630002
C         EIGENVALUES AND EXPANSION COEFFICIENTS - D                    01640000
C                                                                       01650002
      IF(N.GT.M) GO TO 4009                                             01660000
         CALL  DMATR(M,CH,NMAX,AMX,-1)                                  01670002
 4009 MTRX=L+1                                                          01680000
      AINT=AMX(MTRX)                                                    01690000
C                                                                       01700013
         CALL  RDCOEF(CH,M,N,-1,AINT,JW,D,ALAMD,PI,SIGMA,DSUM,FCTR)     01710011
C
      IF( DABS(DSUM).LT.1.0D-50) GO TO 3                                01720000
        IF(INDX.EQ.1)  THEN                                             01730011
         JFN= 2*(JW/2)+1                                                01731018
C+       WRITE(6,999)  CH,M,N,JW,ALAMD,D(1),D(JFN),PI,SIGMA,DSUM,FCTR   01740022
C+999  FORMAT(1H , F7.2,3I5,1P7D14.4 )                                  01750022
        END IF                                                          01760011
C                                                                       01770011
CC            RADIAL FUNCTION AND ITS DERIVATIVE - RJ(M,N,HX),DRJ       01780002
CC            THE RADIAL FUNCTIONS OF FIRST KIND                        01790002
      IFIN=JW                                                           01800000
      IF(MOD(L,2).EQ.0) IFIN=JW+1                                       01810000
      R1A= 0.D0                                                         01820000
      DR1A=0.D0                                                         01830000
      DO 1060 JR=1,IFIN,2                                               01840000
      J=JR-1                                                            01850000
      IF(MOD(L,2).EQ.1) J=JR                                            01860000
      MR=M+J                                                            01870000
      IF(MR.GT.NSJ) GO TO 1060                                          01880000
      IR=IABS(J-L)                                                      01890000
      Q=-1.D0                                                           01900000
      IF(MOD(IR,4).EQ.0) Q=1.D0                                         01910000
      RFAC1= 1.0D0                                                      01920016
      RFAC2= 1.0D0                                                      01930016
        IF(M.GE.1)  THEN                                                01940013
      DO 15 IP=1, M                                                     01950013
      RFAC1= RFAC1*DFLOAT(J+IP)                                         01960013
   15 RFAC2= RFAC2*DFLOAT(J+IP+M)                                       01970013
        END IF                                                          01980013
      RCONV1= D(JR)*RFAC1*DSJ(MR+1,INDX)*RFAC2/FCTR                     01990019
      R1A=R1A+Q*D(JR)*RFAC1*SJ(MR+1,INDX)*RFAC2/FCTR                    02000019
      DR1A=DR1A+Q*RCONV1                                                02010000
      IF(J.GE.IFIN/3.AND.DABS(RCONV1/DR1A).LT.1.0D-14) GO TO 5          02020016
 1060 CONTINUE                                                          02030000
    5 R1=XIA*R1A/PI                                                     02040019
      DR1=XIB*R1+CH*XIA*DR1A/PI                                         02050019
      IF(INDX.EQ.2) GO TO 4                                             02060022
C                                                                       02072015
C          RADIAL FUNCTION OF SECOND KIND IN TERMS OF SPHERICAL NEUMANN 02080002
C                                                                       02090002
      IF(XI.LT.1.005) GO TO 7                                           02100000
      IF(HXX.LT.15.0) GO TO 7                                           02110000
      IF(L.GT.LMTN) GO TO 7                                             02120000
      R2A= 0.D0                                                         02130000
      DR2A=0.D0                                                         02140000
      DO 1070 JR=1,IFIN,2                                               02150000
      J=JR-1                                                            02160000
      IF(MOD(L,2).EQ.1) J=JR                                            02170000
      MR=M+J                                                            02180000
      IF(MR.GE.MAXSN) GO TO 1070                                        02190000
      IR=IABS(J-L)                                                      02200000
      Q=-1.                                                             02210000
      IF(MOD(IR,4).EQ.0) Q=1.                                           02220000
      RFAC3= 1.D0                                                       02230013
      RFAC4= 1.D0                                                       02240013
        IF(M.GE.1)   THEN                                               02250013
      DO 17 IP=1, M                                                     02260014
      RFAC3= RFAC3*DFLOAT(J+IP)                                         02270013
   17 RFAC4= RFAC4*DFLOAT(J+IP+M)                                       02280014
        END IF                                                          02290013
      RCONV2= D(JR)*RFAC3*SN(MR+1,INDX)*RFAC4/FCTR                      02300019
      R2A=R2A+Q*RCONV2                                                  02310000
      DR2A=DR2A+Q*D(JR)*RFAC3*DSN(MR+1,INDX)*RFAC4/FCTR                 02320019
      IF( DABS(RCONV2/R2A).LE.1.0D-08) GO TO 8                          02330000
 1070 CONTINUE                                                          02340000
    8 R2N=XIA*R2A/PI                                                    02350019
      DR2N=XIB*R2N+CH*XIA*DR2A/PI                                       02360019
      R2N=R2N/SCLSN                                                     02370000
      DR2N=DR2N/SCLSN                                                   02380000
      DR2NC=(1.D0/(CH*XI2)+R2N*DR1)/R1                                  02390000
      CRN= DABS((DR2NC-DR2N)/DR2NC)                                     02400000
      R2M(1)= R2N                                                       02410000
      DR2M(1)=DR2N                                                      02420000
      CR(1)= CRN                                                        02430000
    7 CONTINUE                                                          02440000
C                                                                       02450002
CC        SPHEROIDAL RADIAL FUNCTIONS OF THE SECOND KIND                02460002
CC        IN TERMS OF THE ASSOCIATED LEGENDRE FUNCTIONS                 02470002
      IF(L.LT.LM1) GO TO 1050                                           02480000
      IF(MOD(L,2).EQ.1) GO TO 3000                                      02490000
      JFIN=2*M                                                          02500000
      JINT=2                                                            02510000
      GO TO 3010                                                        02520000
 3000 JFIN=2*M-1                                                        02530000
      JINT=1                                                            02540000
 3010 QSUM=(0.D0,0.D0)                                                  02550000
      DQSUM=(0.D0,0.D0)                                                 02560000
      DO 10 I=1,IFIN,2                                                  02570000
      IA=I-1                                                            02580000
      IF(MOD(L,2).EQ.1)  IA=I                                           02590000
      MA=M+IA                                                           02600000
      IF(MA.GE.MAXPQ) GO TO 11                                          02610000
      QD=D(I)*DQX(MA+1)                                                 02620000
      QSUM=QSUM+D(I)*QX(MA+1)                                           02630000
      DQSUM=DQSUM+QD                                                    02640000
      IF(CDABS(QD/DQSUM).LT.1.0D-12)  GO TO 11                          02650008
   10 CONTINUE                                                          02660000
   11 IF(M.EQ.0) GO TO 3030                                             02670000
      DNRAT(JFIN+2)= 0.D0                                               02680000
      DO 20 J=JINT,JFIN,2                                               02690000
      JA=JINT+JFIN-J                                                    02700000
   20 DNRAT(JA)=-E(-JA+2)/(F(-JA)+G(-JA-2)*DNRAT(JA+2))                 02710000
      DN0=D(1)                                                          02720000
      DO 30 JB=JINT,JFIN,2                                              02730000
      DN(JB)=DN0*DNRAT(JB)                                              02740000
      MJB=M-JB                                                          02750000
      MJC=MJB+1                                                         02760000
      IF(MJB.LT.0) GO TO 2050                                           02770000
      QNA=QX(MJB+1)                                                     02780000
      QNB=QX(MJC+1)                                                     02790000
      GO TO 2060                                                        02800000
 2050 MJD=-MJB                                                          02810000
      IF(MJD.EQ.1) GO TO 2070                                           02820000
      QNA=QN(MJD)                                                       02830000
      QNB=QN(MJD-1)                                                     02840000
      GO TO 2060                                                        02850000
 2070 QNA=QN(1)                                                         02860000
      QNB=QX(1)                                                         02870000
 2060 DQN=(FLOAT(MJB-M+1)*QNB-FLOAT(MJB+1)*XX*QNA)/(XX*XX-1.0)          02880000
      DQN=CI*DQN                                                        02890000
      QSUM=QSUM+DN(JB)*QNA                                              02900000
      DQSUM=DQSUM+DN(JB)*DQN                                            02910000
      DN0=DN(JB)                                                        02920000
   30 CONTINUE                                                          02930000
      IF(MOD(L,2).EQ.1) GO TO 3040                                      02940000
      COEF2=FF(M,N)*FACTMM(M)*DN(JFIN)*C*PI/(DFLOAT(M+M-1)*C**M)        02950012
     +      *FCTR                                                       02960012
      GO TO 3050                                                        02970000
 3040 COEF2=-FF(M,N)*FACTMM(M)*DN(JFIN)*C2*PI/(DFLOAT((M+M-3)*(M+M-1))* 02980012
     +       C**M)*FCTR                                                 02990012
 3050 CONTINUE                                                          03000000
      GO TO 3060                                                        03010000
 3030 IF(MOD(L,2).EQ.1) GO TO 3070                                      03020000
      COEF2=-FF(M,N)*D(1)*C*PI                                          03030000
      GO TO 3060                                                        03040000
 3070 COEF2=-FF(M,N)*D(1)*C2*PI/3.0D0                                   03050000
 3060 CONTINUE                                                          03060000
      IH=H                                                              03070000
      IF(MOD(L,2).EQ.1) GO TO 3080                                      03080000
      KINT=M+M+2                                                        03090000
      KFIN=30+2*(M+IH)                                                  03100000
      GO TO 3090                                                        03110000
 3080 KINT=M+M+1                                                        03120000
      KFIN=31+2*(M+IH)                                                  03130000
 3090 IF(KFIN.GE.138) KFIN=138+MOD(L,2)                                 03140006
 2040 DPRAT(KFIN+2)= 0.D0                                               03150000
      DP0=D(1)                                                          03160000
      IF(M.NE.0)  DP0=DN(JFIN)                                          03170000
      DO 40 KA=KINT,KFIN,2                                              03180000
      KB=KINT+KFIN-KA                                                   03190000
      IF(KB.EQ.KINT) GO TO 2010                                         03200000
      DPRAT(KB)=-E(-KB+2)/(F(-KB)+G(-KB-2)*DPRAT(KB+2))                 03210000
      GO TO 40                                                          03220000
 2010 IF(MOD(L,2).EQ.1) GO TO 2020                                      03230000
      DPRAT(KB)=C2/FLOAT((M+M+1)*(M+M-1))/(F(-KB)+G(-KB-2)*DPRAT(KB+2)) 03240000
      GO TO 40                                                          03250000
 2020 DPRAT(KB)=-C2/FLOAT((M+M-1)*(M+M-3))/(F(-KB)+G(-KB-2)*DPRAT(KB+2))03260000
   40 CONTINUE                                                          03270000
      PSUM=(0.D0,0.D0)                                                  03280000
      DPSUM=(0.D0,0.D0)                                                 03290000
      DO 50 KC=KINT,KFIN,2                                              03300000
      DP(KC)=DP0*DPRAT(KC)                                              03310000
      KD=KC-M-1                                                         03320000
      KE=KD+1-M                                                         03330000
      IF(KE.GT.MAXPQ) GO TO 50                                          03340000
      PSUM=PSUM+DP(KC)*PX(KE)                                           03350000
      DPSUM=DPSUM+DP(KC)*DPX(KE)                                        03360000
      DP0=DP(KC)                                                        03370000
   50 CONTINUE                                                          03380000
          WRITE(6,555)  (KD,DP(KD),KD=KINT,KFIN,2)                      03381027
  555     FORMAT(1H / 8(I5,1PE11.3))                                    03391027
      KFINM=KFIN-M-M                                                    03392027
      IF(KFINM.GE.MAXPQ) GO TO 2030                                     03400000
      IF(KFIN.GE.137) GO TO 2030                                        03410006
      CRIT=CDABS(DP(KFIN)*PX(KFINM)/PSUM)                               03420000
      KFIN=KFIN+4                                                       03430000
      PSCRT=(1.0E-12)*(10.0**(L/10.0))                                  03440000
      IF(CRIT.GT.PSCRT) GO TO 2040                                      03450000
 2030 KFIN=KFIN-4                                                       03460000
      R2Q=(QSUM+PSUM)/COEF2                                             03470000
      DR2Q=(DQSUM+DPSUM)/COEF2                                          03480000
      DR2QC=(1.D0/(CH*XI2)+R2Q*DR1)/R1                                  03490000
      CRQ= DABS((DR2QC-DR2Q)/DR2QC)                                     03500000
      R2M(2)= R2Q                                                       03510000
      DR2M(2)=DR2Q                                                      03520000
      CR(2)= CRQ                                                        03530000
 1050 CONTINUE                                                          03540000
C                                                                       03550002
C            THE RADIAL FUNCTION OF THIRD KIND BY HANISH,S METHOD       03560000
C                                                                       03570002
      IF(H.LT.7.0) GO TO 1080                                           03580000
      IF(L.GT.LMTN) GO TO 1080                                          03590000
      LMT=3*(H+ DABS(ALAMD)/H)                                          03600000
      NUP=LMT+20                                                        03610000
      ISTP=MAXPQ-2                                                      03620000
      IF(NUP.LE.ISTP) NUP=ISTP+10                                       03630000
      IF(NUP.GE.499) NUP=499                                            03640000
      DO 55 I=1,ISTP                                                    03650000
   55 A(I)= 0.D0                                                        03660000
      MP1=M+L+1                                                         03670000
      UP(1)=-3.0*(ALAMD+CH2)/(2.0*CH*DFLOAT(M+1))                       03680000
      DO 60 IR=1,MP1                                                    03690000
      I=IR-M                                                            03700000
      MI=M+I                                                            03710000
      UP(IR+1)=DFLOAT(2*MI+3)*(DFLOAT(MI*(MI+1))-ALAMD-CH2) /(2.0*CH*   03720000
     +  DFLOAT((MI+1)*(MI+M+1)))+DFLOAT(2*MI+3)*DFLOAT(I*MI)/(UP(IR)*   03730002
     +  DFLOAT(2*MI-1)*DFLOAT((MI+1)*(MI+M+1)))                         03740002
   60 CONTINUE                                                          03750000
      A(1)=UP(1)                                                        03760000
      DO 62 I=1,MP1                                                     03770000
   62 A(I+1)=A(I)*UP(I+1)                                               03780000
      DW(NUP+1)=(CH+ALAMD/CH)/DFLOAT(NUP-M)                             03790000
      DO 64 JR=MP1,NUP                                                  03800000
      J=(NUP+MP1)-JR-M                                                  03810000
      MJ=M+J                                                            03820000
      IF(J.LE.0) GO TO 64                                               03830000
      DW(MJ)=2.0*CH*DFLOAT(J*MJ)/DFLOAT(2*MJ-1)/(2.0*CH*DFLOAT((MJ+1)*  03840000
     + (MJ+M+1))*DW(MJ+1)/DFLOAT(2*MJ+3)-(DFLOAT(MJ*(MJ+1))-ALAMD-CH2)) 03850002
   64 CONTINUE                                                          03860000
      DWRAT=UP(MP1)/DW(MP1)                                             03870000
      ISTP1=ISTP+1                                                      03880000
      DO 66 I=MP1,ISTP1                                                 03890000
   66 DW(I)=DWRAT*DW(I)                                                 03900000
      DO 68 K=MP1,ISTP                                                  03910000
      IF( DABS(A(K)).LT.1.0D-50) GO TO 69                               03920000
      IF( DABS(A(K)).GT.1.D+70) GO TO 69                                03930000
   68 A(K+1)=A(K)*DW(K+1)                                               03940000
   69 CONTINUE                                                          03950000
      B=CI**(M+M+1)/(CH*FACT(M))*CDEXP(CI*(CH*XI-DFLOAT(N+1)*PAI2))     03960000
      HSUM=QX(1)                                                        03970000
      DHSUM=DQX(1)                                                      03980000
      DO 70 I=1,ISTP                                                    03990000
      HCONV=A(I)*DQX(I+1)                                               04000000
      HSUM=HSUM+A(I)*QX(I+1)                                            04010000
      DHSUM=DHSUM+HCONV                                                 04020000
      IF(CDABS(HCONV/DHSUM).LE.1.D-18) GO TO 72                         04030000
   70 CONTINUE                                                          04040000
   72 R3=B*HSUM                                                         04050000
      DR3=CI*CH*R3+B*DHSUM                                              04060000
      R2H=(R3-R1)/CI                                                    04070000
      DR2H=(DR3-DR1)/CI                                                 04080000
      DR2HC=(1.D0/(CH*XI2)+R2H*DR1)/R1                                  04090000
      CRH= DABS((DR2HC-DR2H)/DR2HC)                                     04100000
      R2M(3)= R2H                                                       04110000
      DR2M(3)=DR2H                                                      04120000
      CR(3)= CRH                                                        04130000
 1080 CONTINUE                                                          04140000
C                                                                       04150002
C         CHOOSSING THE BEST VALUES.                                    04160000
C                                                                       04170002
      SMLST=DMIN1(CR(1),CR(2),CR(3))                                    04180002
      DO 7000 II=1, 3                                                   04190000
      I=II                                                              04200000
      IF(CR(II).LT.1.D-13) GO TO 7010                                   04210000
      CMETHD=DABS((SMLST-CR(II))/CR(II))                                04220000
      IF(CMETHD.LT.1.D-7) GO TO 7010                                    04230000
 7000 CONTINUE                                                          04240000
      R2= R2M(2)                                                        04250000
      DR2=DR2M(2)                                                       04260000
      IF(L.GT.MTHD) GO TO 7020                                          04270000
      R2= R2M(3)                                                        04280000
      DR2=DR2M(3)                                                       04290000
      GO TO 7020                                                        04300000
 7010 R2= R2M(I)                                                        04310000
      DR2=DR2M(I)                                                       04320000
C+        WRITE(6,777) R2, I                                            04321022
C+777     FORMAT(1H+, 63X, 'R2',1PE15.4, I3 )                           04322022
 7020 CONTINUE                                                          04330000
      GO TO 4                                                           04340000
    3 WRITE(6,100)                                                      04350000
  100 FORMAT(1H0,20X,'*  ROBSWF  -  RETURN WITHOUT CALCULATION  *  COEFF04360000
     1ICIENTS D ARE ALL ZERO' )                                         04370000
    4 RETURN                                                            04380000
      END                                                               04390000
